#if ABS_SURFACE
#define DELTA_EPS		1.0e-4
#define	DELTA_EPS_INV	10000

vec3 GetBernstein2(float u)
{
	float cu = 1.0f-u;
	//Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u;
	return vec3(cu*cu,2*u*cu,u*u);
}

float EvalAlgSurf(float u, float v, float w)
{
	vec3 nn;
	float f1[9];
	float f2[3];
	nn = GetBernstein2(w);
	// 	f1[0][0] = func[ 0]*nn[0]+func[ 1]*nn[1]+func[ 2]*nn[2];
	// 	f1[0][1] = func[ 3]*nn[0]+func[ 4]*nn[1]+func[ 5]*nn[2];
	// 	f1[0][2] = func[ 6]*nn[0]+func[ 7]*nn[1]+func[ 8]*nn[2];
	// 	f1[1][0] = func[ 9]*nn[0]+func[10]*nn[1]+func[11]*nn[2];
	// 	f1[1][1] = func[12]*nn[0]+func[13]*nn[1]+func[14]*nn[2];
	// 	f1[1][2] = func[15]*nn[0]+func[16]*nn[1]+func[17]*nn[2];
	// 	f1[2][0] = func[18]*nn[0]+func[19]*nn[1]+func[20]*nn[2];
	// 	f1[2][1] = func[21]*nn[0]+func[22]*nn[1]+func[23]*nn[2];
	// 	f1[2][2] = func[24]*nn[0]+func[25]*nn[1]+func[26]*nn[2];

	f1[0] = func[0].x*nn[0]+func[0].y*nn[1]+func[0].z*nn[2];
	f1[1] = func[0].w*nn[0]+func[1].x*nn[1]+func[1].y*nn[2];
	f1[2] = func[1].z*nn[0]+func[1].w*nn[1]+func[2].x*nn[2];
	f1[3] = func[2].y*nn[0]+func[2].z*nn[1]+func[2].w*nn[2];
	f1[4] = func[3].x*nn[0]+func[3].y*nn[1]+func[3].z*nn[2];
	f1[5] = func[3].w*nn[0]+func[4].x*nn[1]+func[4].y*nn[2];
	f1[6] = func[4].z*nn[0]+func[4].w*nn[1]+func[5].x*nn[2];
	f1[7] = func[5].y*nn[0]+func[5].z*nn[1]+func[5].w*nn[2];
	f1[8] = func[6].x*nn[0]+func[6].y*nn[1]+func[6].z*nn[2];

	nn = GetBernstein2(v);
	f2[0] = f1[0]*nn[0]+f1[1]*nn[1]+f1[2]*nn[2];
	f2[1] = f1[3]*nn[0]+f1[4]*nn[1]+f1[5]*nn[2];
	f2[2] = f1[6]*nn[0]+f1[7]*nn[1]+f1[8]*nn[2];

	nn = GetBernstein2(u);
	return f2[0]*nn[0]+f2[1]*nn[1]+f2[2]*nn[2];
}
vec3 EvalAlgSurfNabla(float u, float v, float w)
{
	vec3 nabla;
	float f = EvalAlgSurf(u,v,w);
	nabla.x = EvalAlgSurf(u+DELTA_EPS,v,w);
	nabla.x = (nabla.x-f)*DELTA_EPS_INV;
	nabla.y = EvalAlgSurf(u,v+DELTA_EPS,w);
	nabla.y = (nabla.y-f)*DELTA_EPS_INV;
	nabla.z = EvalAlgSurf(u,v,w+DELTA_EPS);
	nabla.z = (nabla.z-f)*DELTA_EPS_INV;	
	return nabla;
}
#else

float EvalAlgSurf(float x, float y, float z)
{
#if _GAS_SPHERE_
	return x*x+y*y+z*z-1.0;							//0 shpere
#endif
#if _GAS_MCMULLENK3_
	return (1+x*x)*(1+y*y)*(1+z*z)+8*x*y*z-2.1;					//4 McMullen K3 model
#endif
#if _GAS_TANGLE_
	float x2=x*x,y2=y*y,z2=z*z; return x2*(x2-5)+y2*(y2-5)+z2*(z2-5)+11.8;				//6 Tangle
#endif
#if _GAS_SCHWARTZ_
	return -(cos(x) + cos(y) + cos(z));  //28 k3DSurf schwartz
#endif
#if _GAS_TRICYLINDER_
	return (x*x+y*y-1)*(x*x+z*z-1)*(y*y+z*z-1)-1;		//29 TriCylinder
#endif
#if _GAS_BLOBSCHWARTZ_
	return   -(x*x+y*y+z*z)+cos(5*x)+cos(5*y)+cos(5*z)-.1;	//30 blobSchwartz
#endif
#if _GAS_DOUBLETORUS_
	float x2=x*x,y2=y*y,z2=z*z; return 0.04-x2*x2*(1-2*x2+x2*x2)+y2*(2*x2*(1-x2)-y2)-z2;  //23 Double torus
#endif
#if _GAS_LemniscateOfGerono_
	return x*x*(x*x-1)+y*y+z*z;														//17 Lemniscate of Gerono
#endif
}
vec3 EvalAlgSurfNabla(float x, float y, float z)
{
	vec3 df;

#if _GAS_SPHERE_
	df.x=2*x; df.y=2*y; df.z=2*z;											//0 shpere
#endif
#if _GAS_MCMULLENK3_
	float x2=x*x,y2=y*y,z2=z*z; df[0]=(1+y2)*(1+z2)*2*x+8*y*z; df[1]=(1+x2)*(1+z2)*2*y+8*x*z;df[2]=(1+x2)*(1+y2)*2*z+8*x*y;		//4 McMullen K3 model
#endif
#if _GAS_TANGLE_
	df.x=x*(4*x*x-10); df.y=y*(4*y*y-10); df.z=z*(4*z*z-10);					//6 Tangle
#endif
#if _GAS_SCHWARTZ_
	df[0]= sin(x); df[1]=sin(y); df[2]=sin(z);  //28 k3DSurf schwartz
#endif
#if _GAS_TRICYLINDER_
	float x2=x*x,y2=y*y,z2=z*z,t=(x2+y2+z2-1),t0=t-z2,t1=t-y2,t2=t-x2; df[0]=2*x*(t0+t1)*t2,df[1]=2*y*(t0+t2)*t1,df[2]=2*z*(t1+t2)*t0;		//29 TriCylinder
#endif
#if _GAS_BLOBSCHWARTZ_
	df[0]=-2*x-5*sin(5*x),df[1]=-2*y-5*sin(5*y),df[2]=-2*z-5*sin(5*z);	//30 blobSchwartz
#endif
#if _GAS_DOUBLETORUS_
	float x2=x*x,y2=y*y; df[0]=4*x*((1-2*x2)*y2-x2*(1-3*x2+2*x2*x2)); df[1]=2*y*(3*x2*(1-x2)-2*y2); df[2]=-2*z; //23 Double torus
#endif
#if _GAS_LemniscateOfGerono_
	df.x=x*(4*x*x-2);df.y=2*y; df.z=2*z;										//17 Lemniscate of Gerono
#endif
	//float x2=x*x,y2=y*y; df.x=4*x*((1-2*x2)*y2-x2*(1-3*x2+2*x2*x2)); df.y=2*y*(3*x2*(1-x2)-2*y2); df.z=-2*z; //23 Double torus
	//df.x=2*x; df.y=2*y; df.z= 3*z*z-2*z;  //27 dingdong surface

	return df;
}


#endif